load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/ut_string.ncl"
load "/public/users/wrf/Codes/WRF_scripts/utils/substring.ncl"
load "/public/users/wrf/Codes/WRF_scripts/utils/WRF_dyy_contributed.ncl"


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;            User Config                   ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; radius to degree
r2d = 45.0/atan(1.0)

if (.not. isvar("configFile")) then
    configFile = "/public/users/wrf/Codes/WRF_scripts/WRF_Post/StationLists/allvalue_stnlist.csv"
    print("configFile is undefined, using default: " + configFile)
else
    print("Using configFile: " + configFile)
end if

lines = readAsciiTable(configFile,1,"string",1)
delim = ","

stnnum    = ndtooned(str_get_field(lines,1,delim))
stnlist    = ndtooned(str_get_field(lines,2,delim))
stnlon    = ndtooned(tofloat(str_get_field(lines,3,delim)))
stnlat    = ndtooned(tofloat(str_get_field(lines,4,delim)))
stndomain = ndtooned(str_get_field(lines,5,delim))
nstns     = dimsizes(stnlist)

levels             = (/10, 30, 50, 70/)
;levels             = (/10, 30, 50, 70, 80, 120/)
nlevels    = dimsizes(levels)

;;;;;;;;;;;;;;;;;;;;;;
; Env and Settings   ;
;;;;;;;;;;;;;;;;;;;;;;


nlevel = dimsizes(levels)

if (.not. isvar("output_prefix")) then
    output_prefix = "~/wrfpost/test/station_allvalue"
    print("output_prefix is undefined, using default = " + output_prefix)
else
    print("Using output_prefix: " + output_prefix)
end if

if (.not. isfilepresent(output_prefix)) then
    system("mkdir " + output_prefix)
end if

; test if input_file is given
if (.not. isvar("input_file")) then
    input_file = "/public/users/wrf/wrfout/2012103120/wrfout_d02_2012-10-31_12:00:00"
    print("input_file is undefined, using default = " + input_file)
else
    print("Using input_file: : " + input_file)
end if

; default output time interval is 15min
if (.not. isvar("output_interval")) then
    output_interval = 0.25 ; default is 15 min
    print("output_interval is undefined, using default = " + output_interval)
else
    print("Using output_interval = " + output_interval)
end if

a = addfile(input_file + ".nc", "r")
filename = systemfunc("basename " + input_file)

; check start_date
if (.not. isvar("start_date")) then
    start_date = substring(filename, 11, 23)
    print("start_date is undefined, using default = " + start_date)
else
    print("Using start_date: " + start_date)
end if

domain = substring(filename, 7, 9)

times = a->Times
hours = WRF_Times2Udunits_c(times, 0)
ntimes = dimsizes(hours)

; UTC -> CST
hours = hours + 8 ; 
timeStr_cst = ut_string(hours, "%Y-%N-%D_%H:%M:%S")
startdate_cst = ut_string(hours(0), "%Y-%N-%D_%H")

wrf_dt = (hours(ntimes-1) - hours(0))/(ntimes-1)
dt = round(output_interval / wrf_dt, 3)

; create time series in desired time resolution
output_hours = fspan(hours(0), hours(ntimes-1), tointeger((hours(ntimes-1) - hours(0))/output_interval + 1))
output_ntime = dimsizes(output_hours)
copy_VarAtts(hours, output_hours)

;;;;;;;;;;;;;;;;;;;;;;;;;
; Loop for each station ;
;;;;;;;;;;;;;;;;;;;;;;;;;

do istn = 0, nstns - 1
    print("Processing " + stndomain(istn))
    if (.not. isStrSubset(filename, stndomain(istn))) then
        print("The domain is not desired, skipping.")
        continue
    end if


    output_filename = output_prefix + "/allvalue_" + stnnum(istn) + "_" + stnlist(istn) +  "_" + domain +  "_" + startdate_cst + ".csv"
    ;;; latlon to XY ;;;
    res = True
    res@returnInt = False   ; return real values
    loc = wrf_user_ll_to_ij(a, stnlon(istn), stnlat(istn), res) - 1

    sxc0 = floattointeger(loc(0))
    syc0 = floattointeger(loc(1))

    syc1 = syc0 + 1
    sxc1 = sxc0 + 1

    print ("X/Y location is: " + loc)

    line_array = new(ntimes(0)*nlevel, "string")
    i_line = 0
    ;    line_array(0) = "time,                   u10,       v10,       ulvl1,     vlvl1,    hlvl1,      ulvl2,     vlvl2,    hlvl2,      alfa_u,    alfa_u2,   alfa_v,    alfa_v2,   u70,       u70_2,     v70,       v70_2"
    ;    i_line = i_line + 1

    ;;; loop for time steps ;;;
    ;    do istep = 0, ntimes(0) - 1


    ;;; Load Terrain Data ;;;
    yy = (/syc0, syc1/)
    xx = (/sxc0, sxc1/)
    kk = (/0, 1, 2/)

    ;;; load variables
    ; read variables from wrfout ncfile
    ter         = wrf_dyy_getvar(a,"ter",0,kk,yy,xx)
    z           = wrf_dyy_getvar(a,"z",-1,kk,yy,xx)
    pres_hpa    = wrf_dyy_getvar(a, "pressure", -1, kk, yy, xx)
    uvmet = wrf_dyy_getvar(a, "uvmet", -1, kk, yy, xx)
    uvmet10 = wrf_dyy_getvar(a, "uvmet10", -1, kk, yy, xx)

    ; interpolate to the site location
    tmp_uv = rm_single_dims(linint2( xx, yy, uvmet, False, loc(0), loc(1), 0))
    tmp_uv10 = rm_single_dims(linint2( xx, yy, uvmet10, False, loc(0), loc(1), 0))
    tmp_ter = rm_single_dims(linint2( xx, yy, ter, False, loc(0), loc(1), 0))
    tmp_z = rm_single_dims(linint2( xx, yy, z, False, loc(0), loc(1), 0))
    tmp_pres = rm_single_dims(linint2( xx, yy, pres_hpa, False, loc(0), loc(1), 0))
    delete(uvmet)
    delete(uvmet10)
    delete(ter)
    delete(z)
    delete(pres_hpa)

    ; get surface pressure
    tmp_surf_pres = tmp_pres(:,0)
    surf_pres = dyy_interp_ts(hours, tmp_surf_pres, output_hours, "cubic")
    delete(tmp_surf_pres)

    ; do the vertical interpolation
    uv_arr = wrf_uv_intrp_vert(tmp_ter, tmp_z, tmp_uv, tmp_uv10, levels)
    delete(tmp_uv)
    delete(tmp_uv10)
    delete(tmp_ter)
    delete(tmp_z)
    delete(tmp_pres)

    ; calculate wind speed and direction
    spd_transpose = transpose(sqrt(uv_arr(0,:,:)^2 + uv_arr(1,:,:)^2))
    dir_transpose = transpose(round(atan2(uv_arr(0,:,:), uv_arr(1,:,:)) * r2d + 180, 3))

    ; interpolate to 15 min
    output_spd = dyy_interp_ts(hours, spd_transpose, output_hours, "cubic")
    output_dir = tointeger(dyy_interp_ts(hours, spd_transpose, output_hours, "direction"))
    delete(spd_transpose)
    delete(dir_transpose)
    delete(uv_arr)
    ;;; END OF WIND EXTRACTION ;;;

    ;;; BEGIN RADIATION EXTRACTION ;;;

    ; interpolate to the site location
    tmp_sswdn = rm_single_dims(linint2( xx, yy, wrf_dyy_getvar(a,"SSWDN",-1, kk, yy, xx), False, loc(0), loc(1), 0))
    sswdn = dyy_interp_ts(hours, tmp_sswdn, output_hours, "linear")
    delete(tmp_sswdn)

    ;;; END OF WIND EXTRACTION ;;;

    ;;; BEGIN PRECIPITATION EXTRACTION ;;;
    ; read precipitaion from ncfile
    rain_exp = wrf_dyy_getvar(a,"RAINNC", -1, kk, yy, xx)
    rain_con = wrf_dyy_getvar(a,"RAINC", -1, kk, yy, xx)
    rain_tot_raw = rain_exp + rain_con
    tmp_rain_tot = rm_single_dims(linint2( xx, yy, rain_tot_raw, False, loc(0), loc(1), 0))
    ; interpolate to 15 min
    rain_tot  = dyy_interp_ts(hours, tmp_rain_tot, output_hours, "linear")
    rain_tend = array_append_record(0.0, rain_tot(1:output_ntime-1) - rain_tot(0:output_ntime-2), 0)
    delete(rain_tot_raw)
    delete(tmp_rain_tot)
    delete(rain_exp)
    delete(rain_con)
    ;;; END OF WIND EXTRACTION ;;;

    ;;; BEGIN TEMPERATURE AND RH EXTRACTION ;;;
    tmp_tc2 = rm_single_dims(linint2( xx, yy, wrf_dyy_getvar(a,"T2",-1,kk,yy,xx)-273.16, False, loc(0), loc(1), 0))
    tmp_rh2 = rm_single_dims(linint2( xx, yy, wrf_dyy_getvar(a,"rh2",-1,kk,yy,xx), False, loc(0), loc(1), 0))

    ; interp to output_interval min resolution
    tc2 = dyy_interp_ts(hours, tmp_tc2, output_hours, "cubic")
    rh2 = dyy_interp_ts(hours, tmp_rh2, output_hours, "linear")
    delete(tmp_tc2)
    delete(tmp_rh2)

    ;;; END OF TEMPERATURE AND RH EXTRACTION ;;;

    ;;; write to txt file ;;;
    header = "ID,  SITE NAME,     TIME"
    ; station number and time
    printArray = stnnum(istn) + "," + stnlist(istn) + "," + ut_string(output_hours, "%Y-%N-%D_%H:%M:%S")
    
    ; wind speed and direction
    do ilvl = 0, nlevel - 1
        header = header + ", " + sprinti("%3i", levels(ilvl)) + "sp, " + sprinti("%3i", levels(ilvl)) + "dir"
        printArray = printArray + ", " + sprintf("%5.2f",output_spd(ilvl,:)) + ", " + sprinti("%3i",output_dir(ilvl,:))
    end do

    ; SSWDN
    header = header + ", sswdn" 
    printArray = printArray + ", " + sprintf("%7.2f",sswdn)

    ; Precipitation
    header = header + ", precipitation, aggregated precipitation" 
    printArray = printArray + ", " + sprintf("%7.4f",rain_tend) + ", " + sprintf("%7.3f", rain_tot)

    ; temp, pres and RH
    header = header + ", temperature, pressure, RH" 
    printArray = printArray + ", " + sprintf("%5.2f",tc2) + ", " + sprintf("%7.3f", surf_pres) + ", " + sprintf("%5.2f", rh2)

    tmp = array_append_record(header, printArray, 0)
    allPrintArray = array_append_record(tmp, "@@@@", 0)
    delete(tmp)

    asciiwrite(output_filename, allPrintArray)
    print("Creating: " + output_filename)

    delete(output_spd)
    delete(output_dir)
    delete(tc2)
    delete(rh2)
    delete(rain_tot)
    delete(rain_tend)
    delete(sswdn)
    delete(printArray)
    delete(allPrintArray)
    delete(surf_pres)
end do
